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This  report  deals  with  the  digital  simulation  of  optical  wavefronts 

Wh\5avi  deqcaded  by  Propagation  through  a turbulent  atmosphere. 

Such  simulated  wavefronts  are  necessary  in  the  computer  simulation 
or  optical  systems  operating  on  atmospherically  degraded  light.  Such 
system  simulation  is  very  desirable  because  of  the  expense  of  building 
prototype  systems  and  testing  them  at  remote  locations  after  construction 
The  prime  motivation  is  computer  testing  of  compensated  imaging  systems 
for  viewing  objects  at  extremely  large  distances  through  the  atmosphere. 

The  computer  generated  wavefronts  to  be  described  here  are  characterized 
by  other  factors.  They  simulate  light  which  has  propagated  vertically 
through  the  atmosphere  and  have  as  a result  an  extremely  long  outer 

and  the  IJfpTnlTSi5  alS°  inClU?6  anplitude  « "ell  ^ phase  fluctuations 
and  the  effect  of  the  cross  correlation  between  them. 

1.  INTRODUCTION 


In  the  past,  two  techniques  have  been  used  to  generate  members 
of  random  ensembles:  the  Fourier  transform  technique  and  the  orthogonal 

polynomial  technique.  In  the  Fourier  transform  technique  the  power 
spectrum  of  the  random  function  is  assumed  known.  An  array  of  uncorrelated 
gaussian  random  numbers  with  zero  mean  and  unit  variance  is  generated 
using  standard  digital  computer  programs.  The  Fourier  transform  of 
the  array  is  computed  and  multiplied  by  the  square  root  of  the  power 
spectrum  and  the  inverse  Fourier  transform  of  the  product  is  computed. 

'he  resu.ting  array  still  has  gaussian  random  variables  but  the  correlation 
function  is  now  that  related  to  the  power  spectrum  as  desired. 

In  the  orthogonal  polynomial  technique  the  random  function  is 
represented  as  a sum  of  orthogonal  polynomials,  with  random  coefficients. 

For  example  one  might  use  Zernike  polynomials  to  represent  a random 
function  in  a round  aperture.  This  scheme  is  frought  with  peril  however, 
unless  the  polynomials  happen  be  statistically  independent.  Otherwise 
t ere  will  be  cross-correlation  between  the  various  polynomial  coefficients. 
A more  desirable  approach  would  use  the  Karhunen-Loeve  technique  (Davenport, 
I958j  m which  a set  of  polynomial s is  generated  which  are  eigenfunctions 
of  the  covariance  function.  The  random  function  is  then  expanded 
as  a series  of  these  eigenfunctions.  Coefficients  of  the  eigenfunctions 
are  then  known  to  be  statistically  independent.  y 


Simultation  of  random  atmospheric  quantities  has  been  performed 
using  the  Fourier  transform  technique  (Hogge,  1973),  (Bradley,  1974)  where 
re  ractive  index  variations  were  generated.  The  Fourier  transform 
technique  has  also  ueen  used  (McGlammery , 1974),  (Brown,  1974)  for  the 
simulation  of  respectively  two  and  one  dimensional  random  phase  fronts. 

Random  wavefronts  have  also  been  generated  using  a sum  of  Zernike 
polynomials  (Noll,  1974)  and  the  Cholesky  decomposition  (Bradley,  1974) 

The  work  herein  also  related  to  the  wavefront  polynomial  fitting  scheme* 
of  Fried  (1965)  and  to  general  work  on  Karhunen-Loeve  series. 
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In  this  report  we  extend  previous  work  in  that  both  log-amplitude 
as  well  as  phase  are  considered  as  well  as  their  cross-correlation. 

We  also  consider  spatial  spectra  of  much  larger  frequency  range  than 
heretofore  used.  The  result  is  an  approach  which  combines  both  the 
polynomial  and  Fourier  transform  technique  approaches  in  an  optimum 
fashion. 

In  this  paper  a random  wavefront  is  represented  as  the  sum  of 
several  terms,  each  pertaining  to  a different  spatial  frequency  range. 

In  each  range  a technique  appropriate  to  the  range  is  used  to  generate 
Doth  phase  and  log-amplitude  fronts.  In  the  highest  spatial  frequency 
range  the  Fourier  transform  technique  is  used,  extended  to  produce 
phase  as  well  as  log-amplitude  fronts.  In  the  lower  spatial  frequency 
* anges  wnere  the  fluctuations  are  less  rapid  a polynomial  scheme  is 
used  which  starts  out  with  a set  of  polynomials  similar  to  Fried's 
but  including  more  terms  and  ends  up  with  the  statistically  independent 
Karhunen-Loeve  polynomials.  These  are  then  summed  with  random  independent 
coefficients  to  provide  a particular  manifestation. 

In  the  balance  of  the  paper  we  will  first  present  the  physical 
picture  considered  along  with  a description  of  the  spectral  separation 
technique.  Then  the  high  frequency  problem  will  be  considered,  in- 
cluding coupled  phase  and  log-amplitude.  In  the  next  section  the 
development  of  the  polynomial  fitting  scheme  used  for  the  lower  spatial 
frequency  ranges  will  be  presented.  This  will  include  the  generation 
of  a set  of  polynomials  and  their  diagonalization  to  form  the  Karhunen- 
Loeve  functions.  Following  that  the  result  of  all  the  techniques 
is  demonstrated  and  a typical  representation  shown.  The  final  section 
contains  discussion,  summary  and  conclusions. 

2.  GENERAL  STATEMENT  OF  PROBLEM  AND 

FREQUENCY  DIVISION 

The  general  physical  picture  which  we  are  considering  is  shown 
in  Fig.  1 where  we  see  a light  beam  propagating  from  a great  height 
down  through  an  atmosphere  with  fluctuating  refractive  index  to  the 
input  aperture  of  an  optical  system.  It  is  assumed  that  the  outer 
scale  L0  and  the  structure  parameter  of  the  turbulent  fluctuations 
vary  with  height. 

We  desire  to  generate  an  ensemble  of  random  phase  and  log-amplitude 
fronts  over  the  input  aperture.  The  phase  and  log-amplitude  fronts 
should  be  normally  distributed  with  the  spatial  autocovariances  and 
crosscovariances  or  equivalently  the  spatial  power  spectra  and  cross 
spectra  appropriate  to  the  situation.  Finally  the  phase  and  log-amplitude 
should  respond  faithfully  to  all  scale  sizes  from  the  largest  to  the 
smallest. 
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Illustration  of  physical  situation  considered. 
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The  spatial  spectra  of  the  phase  and  log-amplitude  auto-covariances 
and  their  cro^s-covariance  are  calculated  in  Appendix  A and  are  shown 
graphically  in  Fig.  2 for  a coherence  length  r0  of  60  cm  wavelength 
of  0.6328  microns  and  ratio  of  source  height  to  receiver  height  of 
lfr.  They  were  calculated  using  standard  techniques  (Tatarski , 1971). 

Examination  of  the  spectra  reveal  one  major  problem,  the  tremendous 
strenqth  associated  with  the  low  spatial  frequencies  corresponding 
to  scale  sizes  much  larger  than  the  aperture.  We  note  that  the  knee 
in  the  phase  spatial  spectrum  curve  occurs  at  = 10~5  rad'rrH  which 

corresponds  to  a spatial  scale  of  6 x 105  meters.  If  we  were  going 
to  use  the  Fourier  transform  technique  without  any  modifications  we 
would  reed  a grid  at  least  10$  meters  on  a side  and  with  elements 
much  smaller  than  the  input  aperture,  say  a millimeter  in  size.  That 
indicates  a square  matrix  array  with  10^  elements  on  a side,  indeed 
not  the  simplest  task.  If  a coarser  grid  were  used  there  would  be 
too  rapid  i change  between  the  lowest  spatial  freouency  component 
and  the  ne>t  higher  one.  With  such  a rapid  change  in  the  spectrum 

the  discrete  Fourier  transform  operation  would  not  faithfully  represent 

the  spatial  behavior  without  aliasing. 

Another  part  of  the  problem  is  that  the  large  scale  fluctuations 

are  much  larger  than  the  aperture  size.  The  large  scale  fluctuations 

will  have  relatively  smooth  variations  over  the  size  of  the  aperture. 
Some  method  must  be  found  for  representing  only  the  portion  of  the 
large  scale  fluctuations  which  are  present. 

In  order  to  solve  the  problem  of  representing  the  large  scale 
fluctuations  we  conceptually  split  the  representation  problem  into 
two  parts:  dividing  the  spatial  spectrum  and  relating  the  portions 
of  the  spatial  spectrum  to  fields  over  the  aperture.  We  first  consider 
the  problem  of  dividing  the  spatial  spectrum  into  sections. 

To  continue  we  split  the  spatial  spectrum  up  into  regions,  each 
region  being  sufficiently  small  so  that  it  can  be  represented  by  a 
reasonable  number  of  discrete  points.  Thus  the  two-dimensional  spatial 
spectrum  would  be 


(1)  F(7)=£  F (7) 

0 n 

where  F0  is  the  contribution  to  the  spatial  spectrum  from  the  highest 
spatial  frequency  range,  F^  is  the  contribution  from  the  next  to  highest 
frequency  range  and  so  on.  The  various  spectral  ranges  are  centered 
about  the  origin  and  have  the  following  frequency  limits. 
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The  fraction  31/32  occurs  because  the  element  (33,33)  was  chosen  to 
represent  the  zero  spatial  frequency.  Thus  each  spectral  reqion 
corresponds  to  a square  spatial  region  32n  meters  on  a side. 

In  the  discrete  representation  used  for  the  spectrum  in  the  highest 
spatial  frequency  range  all  but  the  zero  frequency  spectral  value  are 
identical  to  the  true  spectral  value 


(2)  F#(V,  ) 


F(.  ..  ) 
' x'  y' 


= 2-(k-33) 
= 2-t  ( 33-L  ) 


1 < K,L  < 64 


» • * 
x*  y 


/ 0 


The  value  at  zero  frequency  is  somewhat  arbitrarily  chosen.  This  is 
necessary  because  the  actual  value  for  some  of  the  spectra  is  vastly  dif- 
ferent from  that  of  the  first  harmonic,  so  that  aliasing  would  result  in 
the  corresponding  position  function  if  it  were  used.  To  avoid  this  a 
particular  procedure  is  used.  We  impose  the  requirement  that  a finite 
Fourier  series  representation  of  the  spectrum  represent  the  spectrum 
not  only  at  the  discrete  points  but  in  between  the  points  also.  To 
check  this  requirement  a value  of  the  spectrum  at  zero  frequency  is 
chosen  approximately  equal  to  the  value  at  the  first  harmonic.  The 
spectrum  is  represented  as  a finite  Fourier  series  and  then  evaluated 
at  non-integral  values  of  K and  L.  For  a reasonable  value  of  the  zero 
frequency  component  the  spectrum  varies  smoothly  from  one  point  to  another. 
For  an  improper  choice  there  are  oscillations  or  ringing  between  the 
discrete  points. 


For  the  next,  n=l , region  of  the  spatial  spectrum  a discrete 
representation  is  also  used.  In  this  case  it  will  be  used  as  an  intermediate 
step.  The  discrete  representation  is  formed  from  the  true  spectrum 
with  the  spectrum  for  n=0  subtracted.  In  essence  this  means  that 
the  spectrum  for  n=0  is  evaluated  at  each  of  the  discrete  points  between 
zero  and  the  first  harmonic 
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(3) 


kx  = 2tt  ( K-33)/32 
wy  = 2tt  ( 33-L  )/32)/32 


1 1 K,  L < 64 


a"d  !l!btr5Cte?  fro"lthe  true  spectrum  at  those  points.  This  procedure 
has  the  advantage  that  the  newly  formed  spectrum  goes  to  zero  at  the 

edge  of  the  spectral  region  thus  rendering  it  automatically  band-limited 
The  zero  frequency  spectral  value  in  the  n=l  regions  determined  n a 
manner  identical  to  that  used  in  region  n=0.  determined  in  a 

used  TSereg?Sndn=?.f°r  ^ ^ SpeCtral  regions  iS  identical  t0  that 

The  procedure  for  generating  random  arrays  over  the  one  meter 

f?em^evfr°m  JhVpec^ra  associated  with  all  but  the  highest  spatial 
frequency  spectral  region  involves  a procedure  using  polynomials.  This 

simpler  app^  F’>St  t0  f°™  *1*  » OH.W  a 

With  the  spectra  decomposed  as  indicated  it  is  a simple  matter  to 

arrav^nf^^ff*"  tr?nsform  Pr<Jcedure  to  generate  five  sets  of  random 
arrays  of  different  sizes  and  to  find  the  contribution  to  an  area  one 
meter  square  in  the  center  of  the  larger  arrays.  One  would  need  sSto 
sort  of  interpolation  procedure  to  go  from  the  discrete  arrays  to  a 
much  smaller  grid  a meter  square.  Then  the  sum  of  all  the  contributions 

direc ^approach  W0Uld  represent  the  wavefront.  This  would  be  the 

thp  .nV-K  !?  alternative  and  operationally  faster  procedure  for  finding 

bUt1°nS  ?rr  he  ilquare  meter  aperture  arising  from  the  lower  9 
spectral  ranges.  The  procedure  as  indicated  previously  is  to  expand 

the  contribution  to  the  phase  front  across  the  aperture  in  a series  of 
polynomials  with  random  coefficients.  The  variances  of  the  coefficients 
are  determined  from  the  spectra  for  the  given  range. 

One  might  envision  the  procedure  as  having  the  following  steps. 
Suppose  we  were  to  use  the  Fourier  transform  procedure  to  calculate  a 
large  number  of  random  wavefronts  for  each  of  the  spectral  regions.  The 
wavefront  contribution  corresponding  to  the  highest  frequency  spectral 
range  is  retained  for  use  as  is.  Then  suppose  we  choose  a one  meter 
square  with  64  points  in  the  center  of  all  random  wavefronts  for  each 
of  the  other  spectral  ranges.  We  then  interpolate  from  the  large  square 
grids  to  the  one  meter  square  grid.  The  interpolation  should  be  per- 
formed assuming  that  the  Fourier  series  representation  for  each  random 
manifestation  was  a continuous  one  and  evaluating  it  at  the  642  points 
in  the  central  one  meter  array.  Alternatively  the  sampling  theorem 
could  be  used  for  the  interpolation.  The  values  at  the  64*  points  for 
each  manifestation  would  be  used  to  find  the  coefficients  for  a polynomial 
expansion  of  the  manifestation.  The  coefficients  would  vary  randomly  from 
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one  manifestation  to  another.  Indeed  since  the  wavefronts  vary  randomly 
with  normal  distribution  and  zero  mean  the  polynomial  coefficients  will 
also  vary  randomly  with  zero  mean.  One  would  find  the  variance  of  thP 
coeffK,enU  for  each  spectral  range  by  averaging  the  sJJIre  of  elch 
coefficient  over  the  set  of  manifestations. 

Thus  one  would  end  up  with  a polynomial  series  for  each  spectral  ranoe 
representing  the  front  over  the  central  square  meter  The  coeffin'pnic 
the  polynomials  would  be  normally  distributed  with  zero  mean  and  known 
npnpra t eS * T°  geJer5te  a"other  random  manifestation  one  would  merely 
?rhfmpj  Vew  -et  °+l  coefJ1cients  standard  random  number  generation 

642  the'procedure^is^guicker!'  °f  ra"d°"'  COefficients  ” >«s  than 

m ni-n  pr°ce(,ur®  to  follow  one  further  step  will  be  added,  sliqhtly 
complicating  the  mathematics.  That  is  that  polynomials  will  be  found  that 
will  simultaneously  represent  both  phase  and  log-amplitude.  Also  the 
Fourier  transform  procedure  that  is  used  directly  for  the  highest  spatial 

frequency  portion  of  the  spectrum  must  be  modified  to  simultaneously  qenerat< 
both  phase  and  log-amplitude  fronts.  wiianeousiy  general 

In  the  next  section  the  extension  of  the  basic  Fourier  transform 
procedure  to  include  coupled  phase  and  log-amplitude  will  be  considered 
In  the  section  to  follow  the  polynomial  approach  used  for  ?he  W regions 
of  the  spatial  spectrum  will  be  derived  in  general  terms. 

3.  HIGH  SPATIAL  FREQUENCY  REPRESENTATION 

We  now  describe  the  procedure  used  to  generate  the  contribution  to  the 
random  wavefronts  from  the  highest  spatial  frequency  region.  The  procedure 
used  is  based  on  the  Fourier  transform  approach  but  is  more  general.  The 
approach  is  extended  to  two  random  fronts,  phase  and  log-amplitude, ’which 
are  correlated  The  basic  expressions  derived  here  will  also  be  used  in 
the  next  section  in  the  derivation  of  the  representation  for  the  lower 
spatial  frequency  contributions.  We  assume  that  the  spectra  of  the  phase 
and  log-amplitude  covariance  functions  F<  KXlKy)  and  Ft(<x,Kv)  respectively 
as  well  as  the  spectrum  of  their  cross-covariance,  F K ' V are 
available.  x»  y1  d,e 

Procedure  will  be  to  postulate  a particular  model  and  then  show 
that  it  can  be  made  to  have  the  desired  properties.  To  begin,  for  a single 
manifestation  we  generate  two  square  matrices,  R-j  ( I ,J)  and  R2(  I ,J ) of 
statistically  independent  samples  (Hogge,  1974)  from  a distribution  which 
is  gaussian  with  zero  mean  and  unit  variance.  Stated  mathematically,  this 


(4) 


<Ri(I.J)Rj<K.L)».4l  . , L 


8 


postal^  ?S';LbHraCket.S  denote  tl,e  “lactation  operator.  We  then 
covariances^can^e  ^*3“^,^*“  aBd 


(5a) 


:(K'o'Uo>  * Fa(KVL,o)Rl(K’L>  + Fb(KVUo)R2(K-L> 


(5b)  ^(K.o,L.0)  = Fa(K.o.L.0)R,(K,L)  + Fb(K.o.U  )R,(K,L)  . 


spatial*  frequency  UU^he’UcJaVT'F'  V*"'*  °f  ‘h* 

'•  *»'  depend^pon 'the'particular  ^ 


spatial  frequency  range 

The  desired  expected  values  are 


(6a) 


<:(K,;o*L,o)>  * * 0 


(6b)  <|t(K.0,K.o)|2>  * F,(Kro,Uo) 

(6')  <^(^..u0)i*(^#.ue)> . .;2r,.(S,u#) 


(6d) 


/,KW  2 ' ’o2  F.(K'o*L'o)  • 


Fa- Fb>  and  Fc  we  substitute 


(7a)  ,o"F.  - 


^ - .,,2 


lFal2<lR,l2-lFbl2<lR2l2>HF  |2f  |Fh|2 


(7b) 


(7c) 


-2 


o V<4**  = FaFc*<lR,!2>t!Fbl2'lR2l2^  ■ F/c*  * |Fb|2 


-2r 


of>.  * * lFcl2<|R,!2>*|Fb|2<|R2|2>  = IF  I2  + IF.  I2. 
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In  Eqs.  (7),  use  has  been  made  of  the  statistical  independence  of  the 
matrix  elements  stated  in  Eq.  (4).  Solving  Eqs.  (7)  for  Fa,  F^,  and 
F , we  obtain  (for  the  case  where  Ffl,  Ffc  and  Fc  are  real), 


(8a) 


F -F 

F = M. 


< iry 2F. +f, 

o'*  $ x.  $ ^ 


( 8b) 


'0J%-2F^F< 


(8c) 


F = j.-h 

c 


‘d'VV. 


These  three  spectra  are  displayed  in  Fig.  3 for  the  complete  spectrum. 
Note  that  the  magnitude  of  Fc  is  plotted  since  it  becomes  negative 
at  f % 10". 

Thus  to  generate  one  sample  wavefront  the  procedure  is  to  generate 
the  two  random  matrices,  Ri(K,L)  and  R2(K,L)  and  then  <t(Kx:0,Lic0) 
and  *(Kk0iU0)  using  Fqs.  (5)  and  then  to  take  the  inverse  transform 

using 

-i2*r(I-l)(K-lMJ-l)(L-l)l 
^ N0 

(9a)  *(Ix  ,Jx  ) ■ <Q  LL  ♦(K«0.L«o)  e 


-i2,r(I-1)(K-l)-(J-1)(L-1)l_ 

N 

(9b)  t(Ix0,Jx0)  = « l ll  /(K<0»L<0)  e 


Equations  (5)  and  (8)  completely  describe  the  desired  model  for 
the  generation  of  random  wavefronts.  The  model  requires  the  use 
of  two  arrays  of  random  numbers,  and  three  spectra  which  are  nonlinear 
combinations  of  phase  and  log-amplitude  and  cross  covariance  power 
spectra.  Wavefronts  generated  by  this  model  have  been  demonstrated 
to  have  the  correct  first  and  second  order  statistics. 
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it 


For  the  case  at  hand  where  the  spectra  must  be  divided  into  regions 
is  the  spectra  Fa,  F^,  and  Fc  which  would  be  divided  since  it 
is  those  spectra  from  which  the  random  spectra  and  random  wavefronts 
are  generated.  Thus  the  procedure  would  be  to  take  the  spectra  in 
Fig.  3 which  have  been  formed  from  the  undivided  log-anm] i tude,  phase 
and  cross  spectra  of  Fig.  2 and  apply  the  spectral  divisions,  and 
find  the  best  zero  frequency  value.  Equations  (3)  (8)  and  (o)  then 
give  the  wavefront  manifestations,  applicable  to  any  spectral  region 


Tn  this  section  we  have  considered  the  generalization  of  the 
Fourier  transform  procedure  to  the  case  of  coupled  phase  and  log- 
amplitude.  The  equations  for  generation  of  the  combined  random  fronts 
have  oeen  derived.  The  procedure  is  applicable  to  all  spectral  regions, 
although  it  will  be  used  for  only  the  highest  frequency  region.  The 
procedure  will.be  combined  with  the  polynomial  approach  for\ise  in 
the  lower  spatial  frequency  regions  in  the  next  section. 


4.  LOW  SPATIAL  FREQUENCY  REPRESENTATIONS 


In  this  section  we  consider  in  detail  the  representation  for 
the  lower  frequency  regions  of  the  spatial  spectrum.  The  development 
is  based  on  a polynomial  representation  and  applies  to  both  phase 
and  log-amplitude  and  their  cross  correlation. 


In  the  following  both  the  phase  and  log-amplitude  will  have  polynomial 
representations.  That  is,  each  manifestation  of  both  the  phase  and 
log-amplitude  will  be  represented  over  the  input  aperture  by  the  polynomial 


(10a)  . (r)  = l . (r) 

n n 


(10b)  :(r)  = l ;n.n(f) 


where  the  ,n(r)  are  a finite  set  of  orthonormal  polynomials  defined 
over  the  input  aperture  and  the  coefficients  vn,  *n  are  gaussian 
random  variables  with  zero  mean  and  variance  yet  to  be  determined. 

The  various  random  coefficients  are  in  general  correlated,  cross- 
correlations occuring  both  among  the  various  log-amplitude  coefficients 
and  among  the  various  phase  coefficients  and  also  occuring  between  log- 
amplitude  and  phase  coefficients.  The  specific  form  of  the  individual 
polynomials  for  one  situation  will  be  shown  in  the  next  section. 

Using  Eqs.  (10)  the  contribution  from  the  lower  spatial  frequency 
ranges  then  has  the  form  y 
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N N 

(11a)  ^(r)  = l trV>(?)  + j Z <nOn(r). 


Our  goal  is  to  find  another  representation  in  terms  of  a new 
set  of  orthonormal  polynomials  c-j(n  whose  coefficients  are  also 
gaussian  random  variables  with  zero  mean  Jut  whose  coeffcients  are 
uncorrelated.  The  complex  log-amplitude  will  then  be  expanded  in 
a series  of  the  new  polynomials 


(lib)  «1(rr)  « ) c1r,1(r) 


This  new  polynomial  series  will  then  be  used  for  the  random  wave- 
front  generation.  One  merely  generates  the  random  coefficients, 
cm  and  suns  the  series  at  the  desired  values  of  r.  The  series  with 
uncorrelated  coefficients  is  desirable  because  standard  computer- 
based  random  number  generation  schemes  generally  produce  only  uncorrelated 

random  numbers. 

The  procedure  for  generating  the  polynomials  whose  coefficients 
will  be  uncorrelated  is  to  diagonalize  the  covariance  matrix  of  the 
multivariate  gaussian  distribution  associated  with  the  orthonormal 
polynomials.  If  there  are  N polynomials  in  the  set,  then  there  will 
be  2N  gaussian  random  variables.  N are  for  the  coefficients  of  the 
log-amplitude  representation  and  N are  for  the  phase.  Thus  the  joint 
probability  density  for  the  2N  random  variables  is  (Davenport,  1958) 


(12) 


p(*l*  W *•*  ?n)  = 

(2*)  Q exp  - j ZZQ<* + ^ ^m^nJ-rntn  + ^^n^m^m 


where  the  matrix  Q = C'1  is  the  inverse  of  the  covariance  matrix, 
including  phase  autocovariances,  log-amplitude  autocovariances,  and 
phase-log-amplitude  cross-covariances.  The  covariance  matrix  is 
defined  by 


(13a)  " <l>n^m> 


(13b) 


Cf  0 
'■nym 


= <t.  t,  > 

n m 
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(13c)  C*n*m  = <enV 


The  calculation  of  the  covariance  matrix  elements  will  be  discussed 
later  in  this  section. 

The  procedure  for  findinq  the  new  random  indeDendent  coeffi- 
cients is  well  known:  transform  to  a new  linear  combination  of  coeffi- 
cients which  diagonalizes  the  quadratic  form  in  the  exponent  of  the 
joint  probability  density  function.  To  illustrate  this  we  use  a 
matrix  notation.  The  quadratic  form,  Q,  in  the  exponent  of  Eq.  (12) 
is  written  in  vector  notation 


(14a)  2Q  = vT$  v“ 


where  the  vector,  v,  is 


04b)  »*  ...  VV2.  •••  1N) 


and  the  matrix  ^ is  the  inverse  of  the  covariance  matrix  as  indicated. 
To  diagonalize  the  quadratic  form  find  vectors,  w,  satisfying  the 
matrix  eigenvalue  equation 


(15)  H w - i B 
n n n 


Then  the  diagonalizing  transformation  is  the  modal  matrix,  FT,  whose 
columns  are  the  orthonormal  eigenvectors,  Vf  . 


(16)  FT1  0 M = a 


A is  the  diagonal  matrix  containing  the  eigenvalues.  Then  the  new 
random  coefficients  cn  are  given  by 


(17)  v = M c 
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where 


(18)  c = (c-| , • • • C£^) . 


The  cn  are  the  desired  random  independent  variables.  They  are  gaussian 
with  zero  mean  because  the  and  s.n  are  gaussian  with  zero  mean.  The 

cn  also  have  variances  x~'  given  by  the  diagonalized  inverse  covariance 
matrix. 

/ 

The  complex  log-amplitude  in  Eq.  (11)  can  be  rewritten  in  terms  of 
the  new  coefficients.  Writing  Eq.  (11)  in  terms  of  a function  vector, 


(19)  <Kr)  = ( c i . *•*  j'4 1 • jv'2  *•*  JC'n ) 


we  have 


(20)  •*1(r)  = “T(r)v 

(21)  = ~T(r)  H C* 

(22)  = rT(T)C  = l C nf.n(T) 
where 

(23)  c(r)  = (C<| (r) , ^,(r),  C2N(iO)  - MT7 


Equation  (22)  is  the  desired  result.  The  first  order  complex 

log-amp^i tude,  i],  is  represented  as  a sum  of  orthonormal  polynomials, 

the  f (r),  with  random  uncorrelated  coefficients,  the  c . 

" n 

We  note  some  interesting  features  in  this  solution.  The 
polynomials  f.n(r)  are  in  general  complex,  since  the  diagonalizing 
transformation  mixes  log-amplitude  and  phase.  The  solution  is  also 
identical  to  a solution  of  the  Karhunen-Lofeve  problem.  This  is 
demonstrated  in  Appendix  B. 
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We  gow  proceed  to  derive  a qeneral  exnres<ien  fn-  fh« 
matrix  C from  which  its  inverse  5 -1  / *he  £°vana"ce 

ration  involves  a several  s?ep  procedure Basic*  der1- 

to  conceptually  generating  a ^efron?  which i c^lJsV large TJZ 
JS i ng  the  Fourier  transform  technique  Thp  ^tatict-trai  _ _ ?. 

of  this  large  wavefront  are  thus  known  from  the  proceeding  section 

nn?0tpmhntne  th?J  We  gei?erate  a?  ensemble  of  wavefronts  and  inter-’ 
po  ate  between  the  matrix  of  points  using  the  samolina  thpnrpm 

representation.  Imagine  that  we  then  take  the  centO?  oOesouare 
meter  of  each  ensemble  member  and  expand  that  in  terms  of  the  nnlv 
noma;  representation  of  Bp  (10).  He  can  theS  det™X  each"  of' 

the  covariances^  **>  * le  m*mber  and  further  can  ^ermine 

r Ulr.  of  the  covariance  matrix  In  Bp  (10) 

a(*e  ^e^1ned  over  a one  square  meter  aperture  Usinn 
the  orthonormal ity  assumed  for  the  polynomials  ' *'  ’ ^ 


>/'n(r)  we  can  write 


(24) 


i(r) 


>n(r)  w(r)dr 


where  w(r)  is  the  aperture  function, 

(25)  w(f)  =1  j - 1/2  < x < 1/2 

| - V2  < y < 1/2  . 

0 otherwise  , 


in  Eqs.  (7)  We  thpn  talp  thQ  and  (Kk0,Lk-0)  as  indicated 

over  the  large  array?  lni,erSe  transfo™  to  obtain  the  phase 


(26) 


j(MxolNxo)  1 


((K-1)(H-1)+(L-1)(N-1) 


_ 2 

~ K 

0 


IIt(Kr0.U0)e 


theorem  MSM  WrUe  3 San,p'1ng 
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(27)  ^ ( x ,y ) = l l ^(Mx  ,Nx JsincTi 

M N 00 


, x 
L 0 


- (M-33) 


sinc^ 


f-  - (33-N) 

o 


where  sinc(u)  sin(u)/u.  Combining  Eqs.  (5),  (24),  (26)  and  (27)  gives 

(2R)  'M  ' [Fa<l<VUo>Rl<KVL‘o>  * Fb(kVUo)R2<KVUo)] 

VkVL'o> 


where 


(29a)  Un(K>o,L.0)  = O H Vn(Hx0,Nx0)e  0 


jjp  ((K-1)(M-1)*(L-1)(N-1)) 


and 


(29b)  Vn(MxQ,Nxo) 


n 


dr  w(r)c-  (r)sincT 


fr  1 

, 7-  - (M-33) 

lL  o 

sinciT 

f-  - ( 33-N  )1 1 

xo  J! 

A similar  expression  can  be  derived  for  the  log-amplitude  coefficients. 


(290  • I pc(K.0, l'0)R,(K.0.L.0)tFb(K,0,l.0)R2(K,0,U0)]Un(K.0.LK0) 


The  covariance  matrix  elements  then  follow  from  Eqs.  (13),  (28), 
and  (29). 


(30a)  C 


❖ $ 

n*n 


= £ £ <iFa(KVUoll  +lFb(KVUol  > uJ,<kVUo> 


VkV'o> 


"2r 


= l l OF  (Kr  ,U JU*(Kk  ,LOU  (Kk  ,U„) 
o o o m c o n o o 
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(30b)  CVm  =<V;>  * | T (|Fc(K,0.U0)|2+|Fb(K,0,U0  ); 


Um(^ntL<„|Un(kK:  fU  ) 

in'  o’  o’  n'  0 0 


(30c) 


'M. 


n"m 


= <£„££*> 

n m 


H {Fa<KVL*o 


)F;(K<  ,L<  ) ♦ |FJK< 


(S’LS>Ui5(KVL^<KVUo> 


Equations  (35)  together  with  Eqs.  (33)  constitute  the  formal  solution 
for  the  covariance  matrix.  In  the  next  section  a particular  set  of 
polynomials  will  be  chosen.  To  go  with  these  the  particular  associated 
Vn(Mx0,Nx0)  and  Un(KK0,U0)  are- evaluated  in  Appendix  C. 

To  summarize,  we  have  considered  in  formal  terms  a method  for  generating 
the  contribution  to  the  random  wavefront  from  the  lower  spatial  frequency 
portions  of  the  spectrum.  A polynomial  representation  is  derived 
in  which  the  complex  log-amplitude  is  expressed  in  terms  of  a series 
of  complex  orthonormal  polynomials  with  random  uncorrelated  coefficients. 

The  variances  of  the  coefficients  are  related  to  the  covariance  matrix 
defined  in  terms  of  the  cross  correlation  between  the  random  expansion 
coefficients.  General  expressions  for  the  covariance  matrix  elements 
are  derived. 
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5.  RESULTS 


We  now  proceed  to  give  the  specific  details  for  the  random  wavefront 
simulation.  First  the  spectrum  is  partitioned  and  the  low  frequency 
problem  is  treated.  Then  after  choosing  a group  of  ten  polynomials 
orthonormal  over  a source  aperture  of  unit  width,  the  resulting  coveriance 
matrices,  eigenvectors  and  eigenvalues  are  listed.  The  high  frequency 
problem  which  has  been  discussed  previously  is  mentioned  only  in 
regards  to  the  superposition  of  the  two  results.  Finally  a typical 
ensemble  member,  a single  wave-front  manifestation,  is  .displayed. 


The  polynomials  employed  in  the  expansion  of  the  wavefront  are 
listed  in  Eqs.  (31) 

Group  I 

Ola) 

(31b) 

(31c) 

Group  II 

(31d)  <l>5(x,y)  = 12  xy 

Group  III 

( 31  e ) ^(x.y)  = 2/3  x 

(31f)  *6(x,y)  ■ 12J1?  (xy2  - jy  x) 

(31g)  *8(x,y)  = 20j7(x3-  ^ x ) 


^(x.y)  = 1 

^3(x,y)  = 6J|'  (x2+y2-  ) 

*4(x,y)  = 6 J|  (x2-y2) 
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Group 

IV 

( 31  h ) 

*2(x»y) 

(31  i) 

v7(*.y) 

(31  j ) 

'9(*,y) 

As 

indicated 

12/l5  (yx2  - jjy) 


20JT  (y3  - y) 


unit  square  as  described  by  Eq.  (25).  urcnonormai  over  a 

The  grouping  of  the  polynomials  in  Eqs.  (31)  was  done  to  indicate 
the  sunilanties  in  symmetry.  The  first  set  C0(r),  o?(r)  and  ipa(r) 
employ  the  same  combinations  used  by  Fried  (1965).  T^ey  are  even 

hh  -ThG  S^°nd  Se- 1 *5(0  is  odd  in  both  * and  y.  The 
rd  set  is  odd  in  x and  even  in  y.  The  polynomials  in  the  fourth 

set  are  identical  with  tnose  in  the  third  set  except  for  a ninety 
degree  rotation.  As  one  might  expect  the  statistical  properties 

of  the  last  two  sets  will  be  identical  because  the  correlation  functions 
and  power  spectra  have  central  symmetry.  unctions 

The  next  step  after  the  choice  of  polynomials  is  the  partitioning 
o the  spectra.  We  used  a trial  and  error  procedure  to  choose  new 

reqio^Fo?*?  ’ Fb’  and,Fc  in  the  highest  spatial  frequency 

region  For  Fg  the  zero  frequency  value  chosen  was  1.75  times  the 

value  for  the  first  harmonic.  For  Fb  and  Fc  the  zero-frequencv  values 

t0  thG  ralUGS  °f  the  '^Pective  first  haSnics"  p\otl 

/ t5Ldrtfri'!cC-  5-eCIra*uSed  t0r  regions  1 and  2 are  shown  in  Figs. 

* circles  indicate  the  values  at  the  discrete  points  while 

the  solid  curve  comes  from  the  interpolated  values.  The  last  curve 
is  an  example  of  a poorly  chosen  zero  frequency  value. 

The  next  step  is  the  generation  of  the  covariance  matrix.  There 
is  a simplification  that  arises  in  the  covariance  matrix  because 
of  the  symmetry  properties  of  the  particular  set  of  polynomials  which 
can  be  most  easily  demonstrated  by  choosing  a particular  order  for 
the  polynomials.  Thus  for  the  function  vector^ r)  in  Eq.  (19)  we 
choose  a complex  twenty  element  vector.  M v ' 
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Fig.  4b,  Linear  plots  of  spectra  Fa,  Fb,  and  Fc  for  spectral  regions 
n=l  and  2.  The  negative  of  Fb  and  Fc  are  plotted. 
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Linear  plots  of  spectra  Fa,  Fb,  and  Fc  for  spectral 
n-1  and  2.  The  negative  of  Fb  and  F£  are  plotted. 
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Fig.  4d.  Linear  plots  of  spectra  Fa,  Fb,  and  Fc  for  spectral  regions 
n=l  and  2.  The  negative  of  F^  and  Fc  are  plotted. 
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k. 


(32) 


^(r)  (jC0(r),j,3(r) J^4(f),jc5(r),jc1(r)lj;.6(r),j^8(f-), 

r ) . 7 ( r ) • j.gfr)  ,i,o(r) 3(r)  ^ (r ) ,<;,.(  r) , 

.l(r),.6{r)„8(r)f,'2(F)l,7(F)l,9(r)). 


The  imaginary  part  is  given  first  for  later 


the  results. 


convenience  in  displaying 


cinnl^e  SiI"pl!ficati0"  is  tha!'  the  covariance  matrix  divides  into 
simpler  submatrices  as  shown  in  Eg.  (33). 


(33)  C = 


6x6 
Group  I 

0 

2 x 2 
II 

0 

6x6 
Group  I I I 

6x6 
Group  IV 

Eq  "(IS)  Mis?sefn™h  Cr?[’a?Ce  ?trix  0f  the  form  shown  *" 
q.  (33J  exists  for  each  of  the  low  frequency  spectral  regions. 

The  covariance  matrices,  eigenvalues  and  eigenvectors  for 
spectral  regions  n=l  and  n=2  are  shown  in  Table  I.  The  submatrices 
are  divided  into  sections.  The  upper  left  and  lower  right  quadrants 
and  1?9"amp1ltude  covariances  respectively.  The  other 
two  quadrants  contain  phase- log-amplitude  cross-covariances. 

The  eigenvalues  are  listed  in  the  next  row  in  decreasing  order 

teres tin^tQ1  not^t^t the • corresP°ndi n9  eigenvector.  It  is  in- 
teresting to  note  that  the  eigenvectors  with  one  excepetion  all 

contain  one  element  very  near  unity  and  with  the  other  elements 
smaller  by  at  least  several  orders  of  magnitude.  This  indicates 
that  for  these  eigenvalues:  first  there  is  very  little  cross - 
coupling  between  phase  and  log-amplitude,  and  second  that  the  poly- 
nomials chosen  are  indeed  quite  reasonable  representations  of  the 
actual  eigenfunctions.  The  one  exception  is  a mixing  between  the 
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consistant  and  spherical  terms  for  n=l.  The  result  is  two  in- 
tensity patterns  neither  of  which  is  constant.  It  seems  reasonable 
that  a constant  intensity  would  not  be  one  of  the  "natural"  patterns 
for  atmospherically  degraded  light. 

The  adjusted  spectra  Fa,  Fb,  and  Fc  for  the  region  n=0  are 
shown  in  Figs.  5.  The  two-dimensional  integral  of  these  are  ob- 
tained and  combined  to  form  the  variances  of  phase  and  log-amplitude. 

These  variances  are  shown  in  Table  2 along  with  the  sums  of  the 
eigenvalues  for  regions  n=l  and  2.  It  is  interesting  that  the  phase 
variances  for  other  than  the  constant  term  are  larger  than  the  phase 
variance  in  region  n=0.  The  log-amplitude  contribution  of  region 
n=0  is  much  larger  than  that  from  region  n=l,  thus  indicating  the 
high  pass'1  property  of  the  atmosphere  to  amplitude  fluctuations. 

Figures  6 are  contour  plots  made  from  a typical  manifestation 
of  a degraded  wavefront.  There  we  see  the  contributions  from  the 
regions  n=2,l  and  0.  The  contour  lines  are  drawn  every  half  wave- 
length for  phase  and  every  half  radian  for  log  amplitude.  As 
expected  the  contours  are  more  closely  spaced  in  the  high  frequency 
regions.  For  lower  spatial  frequencies  than  n=2  the  contribution  is 
merely  a constant  and  is  not  shown. 

The  results  of  Fable  II  tend  to  indicate  that  the  Fourier  trans- 
form technique  and  tht  polynomial  approach  are  both  necessary  and 
important  in  the  complete  description  of  light  that  has  come  down 
through  the  atmosphere.  One  might  ask  if  the  polynomial  approach 
could  be  used  to  alone  simulate  atmospherically  degraded  wavefronts. 

To  answer  that  question  we  consider  some  further  results  of  a 
slightly  different  nature.  The  Karhunen-Loeve  integral  equation 
was  solved  numerically  using  the  correlation  function  derived  from 
the  phase  spectrum  of  Fig.  2.  A grid  of  eight  by  eight  points  was 
used.  The  eigenvalues  are  shown  in  Appendix  D along  with  the  designations 
of  a few  low  order  eigenfunctions.  Reference  to  that  appendix  shows 
that  the  eigenvalues  are  very  closely  spaced  in  value  indicating 
that  it  would  be  difficult  to  reject  any  of  the  eigenfunctions.  Thus 
all  sixty-four  eigenfunctions  would  be  necessary  for  the  coarse  grained 
eight  by  eight  simulation.  Even  if  all  the  eigenfunctions  were  used 
they  would  have  to  te  summed  at  each  point  with  random  coefficients 
to  generate  the  simulated  wavefront.  If  all  those  random  coefficients 
are  going  to  be  necessary  it  seems  much  more  economical  of  time  to 
use  the  Fourier  transform  technique  to  generate  the  wavefronts  directly 
than  to  have  to  sum  over  each  eigenfunction  for  every  point  desired. 

Thus  for  the  high  spatial  frequency  region  the  Fourier  transform 
approach  seems  much  more  economical  than  the  eigenfunction  approach. 

On  the  other  hand  the  polynomial  approach  requires  many  fewer  terms 
in  the  lower  spatial  frequency  region.  The  general  conclusion  then 
is  that  the  combination  is  the  appropriate  one  where  both  high  and 
low  frequencies  are  involved  all  together. 
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TABLE  II 


n 

Variance 

Phase  Log-amplitude 

0 

.2220  . 3017xl0"3 

1 

.5543xl02  .9330x10"® 

2 

.2068x10®  .2788x10"® 

*?ne  °Jher  •t?m  might  be  mentioned  in  passing.  In  the  aooliratinn 
TMs  r4Slre?^Ttaherecov“XdnS 

£•?■=.  Wa  z sss  as  M/ar1- 

on  this  point  ‘ otv,ous.  More  work  could  be  done 


after  the  new  spectra  have  been  computed  for  reqlon  n=0 
1 the  eigenfunctions  and  eigenvalues  have  been  obtained  for  reoioiis 
4 K'  can  actually  generate  the  random  wavefronfs  This  ?s 
done  using  Egs.  (5)  and  (9)  for  the  highest  spectral  region  an5  Eo 

n h»OrnthV0W<T  "S''0"5'  Th*  PrO«du«  iS^ 

in  the  computer  algorithm  shown  in  Fig.  7.  5iratea 

To  summarize,  in  this  section  the  complete  wavefront  cimniat-inn 
scheme  was  illustrated  for  the  spectra  Zer  conMderaJionl  i se? 

° p°lJynomia^  was  chosen  for  the  low  frequency  spectral  reaions 
and  the  covariance  matrix  diagonalized  to  give  thrLrhuneSolve 
polynomials  appropriate  for  uncorrelated  random  coefficients  A 

complete  wave-front  was  illustrated  in  terms  of  the  contributions 
from  the  various  spectral  regions.  tuniriounons 


6.  SUMMARY 

H0nJHJUmm?^iZ?VWe  ha.e  developed  a scheme  for  simulating  randomly 

the  effprtPifCtLieamS  1ncludlng  both  phase  and  log-amplitude  and  * 
Ha6  °f  bheir  cr°ss-correlation.  The  scheme  works  for  optical 

i!rnada*il0ns*hh0'e  S*a  e 1engths  range  from  much  smaller  to  much 
larger  than  the  input  aperture  size.  An  extension  of  the  well-known 
Fourier  transform  method  is  used  for  the  small  scale  fluctuations 
and  a polynomial  approach  is  used  for  the  larger  scale  fluctuations. 
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In  the  report  Section  2 contains  a discussion  of  the  physical 
situation.  This  is  characterized  by  a large  range  of  scale  sizes 
for  the  fluctuations;  a difficult  problem  for  the  Fourier  transform 
procedure  alone  because  a prohibitively  large  array  would  be  required. 
The  spatial  spectrum  is  split  into  regions,  the  contribution  to  the 
degraded  optical  beam  from  each  spectral  region  to  be  computed 
separately.  Treated  in  Section  3 is  the  highest  spatial  frequency 
region  in  which  the  Fourier  transform  approach  is  used.  In  this  section 
the  mathematical  formalism  for  using  the  Fourier  transform  approach 
to  generate  both  phase  and  log-amplitude  fluctuations  including  their 
correlation  is  derived.  In  Section  4 the  contributions  from  the 
lower  spatial  spatial  frequency  regions  are  considered.  These  are 
generated  using  a polynomial  approach.  The  procedure  for  findinq  the 
set  of  polynomials  which  could  be  used  with  the  uncorrelated  random 
numbers  that  can  be  easily  generated  by  a digital  computer  is  worked 
?ut;u  ^ shown  in  Appendix  B that  the  procedure  is  formally  identical 
to  the  harhunen-Loeve  procedure. 


In  Section  5 the  procedures  developed  in  Sections  2,  3,  and  4 
are  demonstrated  in  the  simulation  of  an  optical  beam  that  has  come 
down  through  the  atmosphere.  The  spatial  spectral  division  is  il- 
lustrated and  a particular  set  of  orthonormal  polynomials  is  trans- 
formed into  the  polynomials  appropriate  for  uncorrelated  random 
coefficients.  The  Fourier  transform  technique  is  used  for  the 
hiqhest  frequency  portion  of  the  spatial  spectrum.  Finally  a 
sample  composite  wavefront  is  illustrated,  showing  the  relative 
contributions  from  the  various  frequency  ranges. 


7.  CONCLUSIONS 

The  problem  considered  initially  in  this  project  was  to  find 
a technique  for  the  computer  simulation  of  radomly  degraded  optical 
beams  including  phase  and  log-amplitude  fluctuations.  In  the  process 
two  approaches  were  to  be  considered,  the  Tourier  transform  approach 
u^ed  previously  by  others  and  the  Karhunen-Loeve  polynomial  approach. 

We  can  conclude  that  a technique  has  been  found.  It  is  to  first  find 
spatial  spectra  that  can  be  combined  to  give  phase  and  log -amplitude; 
second  to  split  these  spectra  into  regions,  each  with  its  own  spectra; 
and  third  to  find  the  contributions  to  the  wave  front  from  each 
spectral  region.  An  extension  of  the  Fourier  transform  technique  is 
used  in  the  highest  spectral  region  while  a polynomial  scheme  is 
used  in  the  lower  frequency  regions.  The  technique  is  one  which  makes 
efficient  use  of  computer  time  because  a polynomial  approach  is  used 
to  simulate  the  effects  of  the  large  scale  fluctuations  which  vary  more 
slowly  over  the  input  aperture  and  the  Fourier  transform  procedure  is 
used  where  it  is  best  suited  in  the  simulation  of  small  scale  fluctuations. 
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During  the  process  of  developing  this  wavefront  simulation  algorithm 
several  other  problems  were  solved.  The  problem  of  using  the  Fourier 
transform  method  to  simulate  both  phase  and  log-amplitude  fluctuations 
,'icluding  their  cross-correlation  was  solved.  The  problem  of  using  the 
Karhunen-Loeve  approach  for  wavefront  simulation  was  considered  and 
extended  to  a complex  variable,  the  complex  log-amplitude,  thus  enabling 
the  application  of  that  approach  to  both  phase  and  loq-amplitude,  A 
problem  of  handling  the  large  range  of  spatial  scale  sizes  was  solved  by 
the  spatial  frequency  division.  Further  the  problem  was  solved  using  an 
efficient  scheme  by  combining  the  Fourier  transform  approach  and  the 
polynomial  approach  for  ranges  where  they  can  be  most  suitably  used, 
fhis  technique  will  be  especially  useful  in  the  computer  examination  of 
active  systems  to  be  used  for  the  compensation  of  atmospheric  degradation 
of  images  (satellite,  etc.)  the  light  for  which  has  propagated  large 
distances  down  through  the  atmosphere. 
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APPENDIX  A 
SPECTRA 


In  this  appendix  we  derive  the  expressions  for  the  spatial  spectra 
of  phase  and  loq-ampl i tude  covariances  and  the  phase-loq  amplitude 
cross  covariance  for  light  that  has  propagated  down  throuoh  the 
atmosphere.  The  expressions  draw  upon  work  in  the  literature 
(Tatarski,  1971)  and  end  up  with  computer  evaluation  of  the  associated 
inteqral  expressions.  In  this  case  the  spectra  are  particularly  sig- 
nificant because  they  apply  to  a light  beam  which  has  propagated  down 
through  the  atmosphere  where  there  is  a much  lamer  low  spatial  fre- 
quency content  than  for  light  propagated  horizontally  near  the  ground, 

The  starting  point  for  the  three  spectra  are  expressions  in  the 
literature  (Tatarski,  1971,  Eqs.  5-27  , 5-33,  45-1-1,15,16,  and  46-22, 
23,25,26), 


? fL  ,-L  j‘  y( Z ‘ -z" ) 

(la)  F 1 ( — T , L ) = * exp  ~2jr—  Fn(kT,z'-z")dz'dz" 

• 0-  0 


I I - i ^ ( ?\  - 7 ' - 7 " ) 

(lb)  F2(-  L)  = -k2  ’ ' exp  — — -2T—  Fn(rT,z'-z")dz'dz" 

•0-0 


(lc)  F = \ [F1  + Re(F2)] 

(ld)  F.  = 1 [F 1 - P.e(F2)] 

(le)  F 4 = I Im  f2 

(lf)  Bj,  )--?,  [ Jo(-r  )F((.T,L).Td.T 
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Oq)  Fn(>-j,z'-z")  = 


“Xronentto?d!  'V & °m  ’ a"d  ^ = !‘x^)  is  the  transverse 

in  ?he  ci,o^>y)tr,.T!sos;?ezrin?4  satin9 

J r?p?-the-PHth’  A'?"{k  is  the  1i9ht  wavelength,  and  *n(,:,z)  is  the  re- 
ractive  index  spatial  spectrum.  We  also  assume  that  th’e  index 

afong^he^ath5  °W  y Varyin9  function  of  man  Position,  (l/2)(z'+z") 

. c AS  Step  in  eva]uating  the  integrals  in  Eqs.  (1)  we  switch 

to  sum  and  difference  coordinates  arranging  the  integration' limits  in 
such  a way  as  to  perform  the  difference coordinate  iStegrSioi 


(2a)  K - z'  - z" 


7-  - + 1 

z - n + 2 


\ 


(2b) 


n = 


Z'  + z" 


Combining  Eqs.  (?) , (la),  and  (Ig),  gives 


I /9 

(3)  F^.L)  = k2  ' d 


2n 

-2n 


dn  e 


iii 

2k 


i k c 

d,z:n^,n^e  1 


+ k‘ 


JL  dn 

r2(L'ri)  eiri 

JL/2  J 

-2(L-n)  i 

drz1n^ir,n^e 


i k r 

z 


thpfclrln3  w ’ ^t^tion,  making  the  substitution  n'=L-C  in 
the  second  integral  and  combining  the  integrals  results  in 
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(4a)  F.(.  ,L)  = 4k^  -idn  d»  sine 

1 1 Jo  J-'  1 


( ' (k  »n)  + * ( ► , L — T) )} 
nv  ’ ' n 

In  Eqs.  (4a),  by  definition,  sinc(x)=sin(x)/x.  Performing  similar 
operations  on  the  F^  integral  gives 


2 fL/r2  f 

(4b)  F„(>t,L)  = -4kz  ndn  d>  ,sinc(2nr  ) 

2 ' J0  J 2 2 

^ - i * j ( L - ■ ) ‘it  Tn) 

'j  Jn(~»n)e  k + :n(.  ,L-r  )e  k j 


We  assume  a von  Karman  type  of  refractive  index  spectrum  where 
both  the  structure  parameter,  Cp  and  the  outer  scale,  Lq  are 
functions  of  height  h given  by  Eqs.  (5) 


,5a,  y.,,)=  0.033 

(5b)  C2(h)  = C2(Ho)  (h/Ho)"4/3 
(5c)  Lo(h)  = Lo(Ho)(h/HQ) 


We  are  considering  downward  propagation  from  height  H-j  to  height  H0 
so  that  the  height  is  related  to  the  path  variable,  n,  by 

h = hl  - (HL-HQ) (n/L)  . 
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Hq  is  the  height  at  the  bottom  of  the  path  so  that 


Combining  Eqs.  (lc) , (Id),  (4a),  (5), 

I-  / 2 

(7a)  Ff  = 2k2  x 0.033  C2(H  ) [ ndn 

Jo 


and  (6)  we  have  for 


and 
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Equations  (7) 
sine  function 
nower  11/6 


can  be  further  simplified  after  examining  the  offset  in 
in  comparison  with  the  width  of  the  term  raised  to  the 


the 


(8) 


2 

' T 
2k 


’ T 


<»' 


*T 


Thus  the  offset  o*  the  sine  function  is  never  significant  and  can 
be  neg  e c tec  . Lgu.it.ons  (7)  then  reduces  to 
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(9c)  F 


= 2k2  x 0.033  C2(H  ) 
n'  o' 


M2 


nd-,  dr  sinc(2nr  ) 

J —'jj  ^ z 
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if  k^nH6^1  *°.have,the  spectral  expressions  expressed  not  in  terms 
of  k and  Cn,  but  in  terms  of  the  coherence  lengh,  rn  to  allow  easier 
comparison  with  the  physical  situation.  9 o to  allow  easier 

FOr  thiS  FW  tw  weve^nKtuITKS  ,S 


(10a) 


Fw  = 


X 0.033  cJ(H  ) [L/?  od't  f*  d, 

' O J -OD 


sinc(2n*2) 

2 


itdth^Vl  ihe^iVf/IUrct,on-  ' ' 2".  always  less  then  ...  the 

jUi:  °A±L*,‘V  1 ' .'’acaese^he  site  of  a typical  refraJtlye 


index  f lurtiiatinS  -o  / ■ , ’ u '•J'H'cai  rciracnve 

ie  have  ’‘2r/,T  1S  a1ways  1ess  lhan  the  ranqe.n.  Thus 


(10b) 


F % 2e-2x  -0.033  C2(H  )hl11/3  fL/ 
w n'  o T 

Jo 


am-'f 


“‘ft  -Hi-)) 


•4/3 


r 

J- 


2n  da-  sine  (2na-  ) 


■ 2xk2  x 0.033  <(h0),;"/3  (3Ho)[,  .(M 


1/3- 


A wave  structure  function 


Ola)  Dw  * 6.88(nr  ) 


5/3 
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corresponds  to  a spectrum 


(Hb) 


6.88  Ul/6)25/3  . ',1/3 
r5/3  2.  x 6-7, 7ir- 

0 5 


Comparing  Eqs.  (10c)  and  (11b)  we  see  that 


.245  -11/3 

5/3  ' 

r 


r 

0 


6.8S1  (11/D2573  

1 (202(6/5):  (l/6)k2x0.033C2(H  )(3H  ) 

n 0 0 


_ 6j88 

? 2/  ' Hn  \]/2 

■2-9'k  (3H0)<(l  -Vi^  j j 


3/5 


,3/5 

1 

! 

i 


In  Eqs.  (9)  substitute 


2k2  x 0.033C2(H  ) . iiSJilrt).?5/!.. 


(6/5)  (l/6)(3Ho)(l 


" V/3 


/H 

M 


5/3 


r 

j ° 


This  is  included  because  it  is  often  simpler  to  work  in  terms  of  r0 
rather  than  k and  Cp.  Equation  (9)  and  (12)  are  the  results  of  this 
writeup.  Figures  1 and  3 in  the  text  are  the  results  of  computer 
evaluations  of  the  expressions  given. 
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APPENDIX  B 


In  this  appendix  we  consider  the  polynomial  representation 
used  for  the  lower  frequency  portions  of  the  spatial  spectrum. 
Specifically  we  derive  the  Karhunen-Loeve  expansion  for  a complex 
function  and  show  that  the  procedure  for  diagonalizing  the  gaussian 
multivariate  distribution  is  equivalent  to  solving  the  Karhunen- 
Loeve  problem. 

To  start  we  postulate  formalism  for  the  Karhunen-Loeve 
expansion  for  a random  complex  function.  We  apply  th  s to  the 
complex  log-amplitude  but  the  development  holds  for  any  complex 
function^  Let  the  random  complex  function  have  real  and  imaginary 
parts  t(r)  and  :(r)  respectively  and  represent  it  by  the  complex 
function  vector 


(Bl)  x(r)  = ( (r)  + j : ( r ) ) 


The  covariance  functions  are  represented  by  a complex  function 
matrix 


(B2) 

where 


C(r' ,r) 


cf  Ar‘  ,r) 
jC  (r\r) 


-JC  ;(r\r) 
C (r\r) 


(B3)  C. . (r.r)  =<{.(?').  :(r)> 


(B3b)  C,.(F,r)  = <t(r)  <(F) 

(B3c)  C^(F.F)  « <i(r' ) «(F)> 

The  spectra  of  C;, (r',F),  C,.(r,,F)  and  C;(r',r)  are  calculated 
in  Appendix  A for  the  case  of  interest.  The  Karhunen-Loeve  integral 
equation  is  formally  stated 
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(B4) 


dr  C(r ‘ ,r)  x(r)  * >x(r) 

J J A 


where  A corresponds  to  the  domain  of  interest  where  e (r ) and  ;(r) 

are  non-zero.  Equation  (B4)  corresponds  to  the  two  coupled  inteqral 
equations  3 


(B5a) 


ff 


ff 


drC  (r',r).(r)  + dr  C .(r',r):(r)  = >^r') 

Jin  ■ * ' 9 


(B5b) 


f f 


A drC  . (r ' ,r) . (r)  + dr  C. . (r 1 ,r) : (r)  = -(r1) 

•JA 


Equations  (B5)  reduce  to  the  standard  Karhunen-Loeve  equations  if 
(r)  and  :(r)  dre  uncorrelated. 

To  solve  Eqs.  (B5)  we  use  a polynomial  series  representation. 


(B6a ) 
(B6b) 


(r)  = I 


on  n 


.‘(r)  = >;  , 

on  n 


(r) 

(r) 


The  zero  subscripts  on 

this  expansion  from  the 

text.  The  orthonormal 

we  will  use  a truncated 

to  the  neglected  terms 

are  indeed  negligible. 

plying  through  by  < (r 

m 


the  * pn  and  :on  are  intended  to  differentiate 
si  mi  la ^expansion  in  Eqs.  (10)  in  the  main 
set  <n(r)  should  be  complete.  In  practice 
series  for  which  the  eigenvalues  corresponding 
are  sufficiently  small  so  that  the  terms  droppe< 
Inserting  solutions  (B6)  into  Eqs.  (B5),  multi- 
) and  integrating  with  respect  to  r'  gives 


(B7a) 


i c 


mn 


on 


' on 
mn 


= ). 


om 


(B7b) 


N 

I C. 

1 y 


mn 


N 


f.  + / C = ) a 

on  ( mn'on  on 


1 
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where 


(B7c) 


r r f r 


dr  dr' 


»r)^n(r)w(r)w(r' ) . 


The  stihscripts  » and  y in  Eq.  (B7c)  stand  for  £ and 
\ ) can  be  combined  into  a single  matrix  eigenvalue 


Equations 

equation, 


(B8)  C"  y = >y 

where  the  matrix  f contains  C r and  r ; and  the  vector  7 
contains  both  k.Q  and  * as  indicated  in  EqsT  (B9) 


(B9a)  C = 


t 

I 


s, 

1 

1 

1 

l 

1 

1 

1 

1 

\ 

c 

”11 

i 

c,  ! 

i.  Z I 
UV  J 

c<, 

CCuv 

cf. 

11 

j S, 
1 

V 

c«.»  ! 
WV  1 

c-,o  / 

UV  / 

^B9t>^  y (f01»  *02*  'nu»  ❖, 


ON*  v01,  '*02 


W 


wi t^setting’the  ^ 'e1'-kn°W"  Kh1‘h  «"* 


(BIO)  |C-  aT|  = 0 
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The  resultant  values  for  the  eigenvalues,  Ap  then  give  the  eigenvectors 
yp  and  the  diagonalizing  matrix  M0  whose  columns  are  the  normalized 
eigenvectors.  The  diagonalizing  transformation  is 

(Bll)  y = MQ  z 


We  now  display  the  connection  with  the  multivariate  distri- 
bution indicated  in  the  text.  Multiplying  Eq.  (B8)  from  the  left 
by  \*i  C"'  gives 

(B12a)  ? y = A-1  y 


where 


(B12b)  Q = C 1 


Equations  (B12)  are  identical  with  Eq.  (15)  in  the  text  if  we  equate 
the  orthonormal  polynomials  .n(r)  and  an(r).  Further  the  eigenvalues 
of  the  $ matrix  in  the  multivariate  distribution  are  merely  the  re- 
ciprocals of  those  of  the  Karhunen-Loeve  problem  and  the  eigenvectors 
yn  and  wn  are  identical. 

This  completes  thp  desired  demonstrations:  the  derivation  of  the 
associated  Karhunen-Loeve  problem  and  the^ demonstrati  on  of  the  equiva 
lence  of  the  multivariate  and  Karhunen-Loeve  approaches. 
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APPENDIX  C 


U,M  n?w  expressions  for  the  functions  Un(K  n,L  n)  and 

V(Mx0,.lx0)  which  were  introduced  in  the  main  test  Pqs.°(29?  in  the 
main  text.  M v ' e 


for  th m ! generality,  we  display  only  the  development 
for  the  phase.  To  begin,  assume  that  the  orthonormal  functions  are 
polynomials  in  x a.  .1  y so  that  in  general  e 


(C-l) 


,(x,y)  = V T 

p q 


q xbyq, 

ynp  q y 


Tf'°n  for  a square  aperture  the  function  V (Mx 
in  tne  form  n'  o’ 


Nxq)  can  be  written 


(C-2) 

V'VV  = 

where 

x /2 

(C-3) 

r 0 

tqOO  ■ 

W2 

’npq  p' 


dx  xq  sinc^ 


xo 


- ( N-33) 


in  t °^i ’ L ' in  Eq‘  can  also  be  wri tten 

in  the  more  detailed  form. 


K-4)  U(K.0.L,0)-n  %q  Tp(^o)T(Lo> 


p q 


where 


(W)  W-'o  j,  l.(Hx0)e’"'<H',)(K'1)  . 

the  func,ions  Uo(K  °,L  »> 
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(C6a) 

Uo<KkVLko>  “ 

T (Kk  )T  (L.  ) 
0 0 0 0 

(C6b) 

U,(Kk„,Kk0)  = 

^3~T,>(Kk  ) 

0 0 1 0 

(C6c) 

U2<KIVLko>  ' 

( C6a ) 

U3(KkVLko>  * 

6j[VE'o)VL'o)*V('.iW 

-f  VK'o>VL'„>] 

(C6e) 

V'VV  = 

! CT„(K.„)T2<K.o)  -t2(k.  o)to(uo)] 

(C6f) 

U6(Kk0,Lk0)  = 

-,2  Tl<K'o,Tl<L'o> 

(C6g) 

U6(,VLko>  = 

12^[T2(K.o)T,(2,.o)  -}IT„(K.0)T1(U0)] 

(C6h) 

WLko>  = 

-iz.TS[t1(k.o)t2(l.o)  -Lt,(k.o)  to(l.o)] 

(C6i ) 

VVV  = 

2q.T[To(K.o)I3(L.0)  -|fT0(K.0)T1(U0)] 

(C6j) 

U9  (Kk0>Lk0)  = 

-20.7  tT3(K.0)T0(L-  0)  - |ffT,(K,0)T0(U0) 

61 


APPENDIX  D 


This  appendix  contains  a listing  of  the  sixty  four  eigenvalues 
of  a solution  for  the  Karhunen-Loeve  integral  equation  assuming  an 
eight  by  eight  grid  for  the  aperture  plane.  The  correlation  function 
corresponded  to  the  phase  function  calculated  in  Appendix  A. 


EIGENVALUES 


. 9999E  0 -(0,0) 
.9220E  -7  -(2,0) 
. 1802E  -7 
.8486E  -8 
.4575E  -8 
. 3201E  -8 
. 2281E  -8 
. 2005E  -8 
. 1 586E  -8 
. 1 461 E -8 
. 1 249E  -8 
. 1071 E -8 
.9478E  -9 


. 1875E  -5  -(0,1) 
•4423E  -7  -(2,1) 
. 1590E  -7 
.8486E  -8 
. 41 62E  -8 
.3015E  -8 
. 2254E  -8 
. 1 958E  -8 
. 1 563E  -8 
. 1461 E -8 
. 1 232E  -8 
. 1065E  -8 
. 87 73E  -9 


. 1875E  -5  -(1,0) 
. 4423E  -7  -(1,2) 
. 1 323E  -7 
. 6235E  -8 
. 4044  E -8 
. 3014E  -8 
.2109E  -8 
.1742E  -8 
. 1 563E  -8 
. 1 344E  -8 
. 1203E  -8 
. 1065E  -8 
.8772E  -9 


. 1 599E 

-6 

-(1 

.1) 

. 2537E 

-7 

-(o 

,3) 

. 8994E 

-8 

.6234E 

-8 

.4043E 

-8 

.2565E 

-8 

. 2034E 

-8 

. 1 712E 

-8 

.1479E 

-8 

. 1 343E 

-8 

. 1199E 

-8 

. 9599E 

-9 

. 8269E 

-9 

. 1 347E  -6  (0,2) 
. 2537E  -7  (3,0) 
. 86 36 E -8 
.4873E  -8 
.3321 E -8 
. 2565E  -8 
.2006E  -8 
. 1586E  -8 
•1474E  -8 
. 1264E  -8 
.1072E  -8 
. 9580E  -9 
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